#Redistricting in R
#Begun 12-11-2008 by Justin Kirkland

# Read in the Data
#################SL Extensions##########################################################
########################################################################################
#Begun 3-26-09
library(foreign)
NCSLdata<-read.dta("NCSLdata.dta")
CASLdata<-read.dta("CalSLset.dta")


###########Separate R and D models###############
#################################################
NCSLdata$turnout<-(NCSLdata$democrats+NCSLdata$republican)/(NCSLdata$NCHouseDem+NCSLdata$NCHouseRep)
CASLdata$turnout<-(CASLdata$assdem+CASLdata$assrep)/(CASLdata$totvote)
marg<-seq(0.5, 1, by=0.001)
NCSLd<-subset(NCSLdata, NCSLdata$party1==0)
NCSLr<-subset(NCSLdata, NCSLdata$party1==1)
model5.NC<-glm(changesl~medinc2+NCHouseDem+NCHouseRep+distance2+percentblack+marginality+medinc2:marginality+turnout, family=binomial(link="logit"), data=NCSLd)
summary(model5.NC)
model5.NCr<-glm(changesl~medinc2+NCHouseDem+NCHouseRep+distance2+percentblack+marginality+medinc2:marginality+turnout, family=binomial(link="logit"), data=NCSLr)
summary(model5.NCr)

se.mfx<-2*sqrt(vcov(model5.NC)[2,2]+vcov(model5.NC)[9,9]*marg^2+2*marg*vcov(model5.NC)[2,9])
se.mfxr<-2*sqrt(vcov(model5.NCr)[2,2]+vcov(model5.NCr)[9,9]*marg^2+2*marg*vcov(model5.NCr)[2,9])
#pdf("NCDemMarg.pdf")
plot(marg, model5.NC$coef[2]+model5.NC$coef[9]*marg, type="l", lty=1, lwd=4, col="black", ylim=c(-0.05, 0.05), ylab="Effect of Income on Redistricting", main="Varying Effect of Income as Democratic Districts Become Safer", sub="North Carolina Model", xlab="Incumbent Proportion of Vote")
lines(marg, model5.NC$coef[2]+model5.NC$coef[9]*marg+se.mfx, lty=2, lwd=2)
lines(marg, model5.NC$coef[2]+model5.NC$coef[9]*marg-se.mfx, lty=2, lwd=2)
abline(0,0, lty=2, lwd=0.5)
#dev.off()
#pdf("NCRepMarg.pdf")
plot(marg, model5.NCr$coef[2]+model5.NCr$coef[9]*marg, type="l", lty=1, lwd=4, col="black", ylim=c(-0.05, 0.05) , ylab="Effect of Income on Redistricting", main="Varying Effect of Income as Republican Districts Become Safer", sub="North Carolina Model", xlab="Incumbent Proportion of Vote")
lines(marg, model5.NCr$coef[2]+model5.NCr$coef[9]*marg+se.mfxr, lty=2, lwd=2, col="black")
lines(marg, model5.NCr$coef[2]+model5.NCr$coef[9]*marg-se.mfxr, lty=2, lwd=2, col="black")
abline(0,0, lty=2, lwd=0.5)
#dev.off()

CASLd<-subset(CASLdata, CASLdata$party1==0)
CASLr<-subset(CASLdata, CASLdata$party1==1)
model5.CA<-glm(slchange~med2+assdem+assrep+near_dist+perAA+marginality+med2:marginality+turnout, family=binomial(link="logit"), data=CASLd)
summary(model5.CA)
model5.CAr<-glm(slchange~med2+assdem+assrep+near_dist+perAA+marginality+med2:marginality+turnout, family=binomial(link="logit"), data=CASLr)
summary(model5.CAr)

se.mfx<-2*sqrt(vcov(model5.CA)[2,2]+vcov(model5.CA)[9,9]*marg^2+2*marg*vcov(model5.CA)[2,9])
se.mfxr<-2*sqrt(vcov(model5.CAr)[2,2]+vcov(model5.CAr)[9,9]*marg^2+2*marg*vcov(model5.CAr)[2,9])
#pdf("CADemMarg.pdf")
plot(marg, model5.CA$coef[2]+model5.CA$coef[9]*marg, type="l", lty=1, lwd=4, col="black", ylim=c(-0.05, 0.05), ylab="Effect of Income on Redistricting", main="Varying Effect of Income as Democratic Districts Become Safer", sub="California Model", xlab="Incumbent Proportion of Vote")
lines(marg, model5.CA$coef[2]+model5.CA$coef[9]*marg+se.mfx, lty=2, lwd=2)
lines(marg, model5.CA$coef[2]+model5.CA$coef[9]*marg-se.mfx, lty=2, lwd=2)
abline(0,0, lty=2, lwd=0.5)
#dev.off()
#pdf("CARepMarg.pdf")
plot(marg, model5.CAr$coef[2]+model5.CAr$coef[9]*marg, type="l", lty=1, lwd=4, col="black", ylim=c(-0.05, 0.05), ylab="Effect of Income on Redistricting", main="Varying Effect of Income as Republican Districts Become Safer", sub="California Model", xlab="Incumbent Proportion of Vote")
lines(marg, model5.CAr$coef[2]+model5.CAr$coef[9]*marg+se.mfxr, lty=2, lwd=2, col="black")
lines(marg, model5.CAr$coef[2]+model5.CAr$coef[9]*marg-se.mfxr, lty=2, lwd=2, col="black")
abline(0,0, lty=2, lwd=0.5)
#dev.off()

##########PredProb
library(mvtnorm)
library(arm)
sims<-500
count.sim<-rmvnorm(sims, model5.NC$coef, vcov(model5.NC))
calc<-400
med<-seq(0, 200, length=calc)
pe<-numeric(calc)
lb<-numeric(calc)
ub<-numeric(calc)
pe2<-numeric(calc)
lb2<-numeric(calc)
ub2<-numeric(calc)
for(i in 1:calc){
	pc1<-numeric(sims)
	pc2<-numeric(sims)
	for(j in 1:sims){
	pc1[j]<-count.sim[j,1]+count.sim[j,2]*med[i]+count.sim[j,3]*mean(NCSLdata$NCHouseDem, na.rm=T)+count.sim[j,4]*mean(NCSLdata$NCHouseRep, na.rm=T)+count.sim[j,5]*mean(NCSLdata$distance2, na.rm=T)+count.sim[j,6]*mean(NCSLdata$percentblack, na.rm=T)+count.sim[j,7]*0.5+count.sim[j,8]*mean(NCSLdata$turnout, na.rm=T)+count.sim[j,9]*0.5*med[i]
	pc2[j]<-count.sim[j,1]+count.sim[j,2]*med[i]+count.sim[j,3]*mean(NCSLdata$NCHouseDem, na.rm=T)+count.sim[j,4]*mean(NCSLdata$NCHouseRep, na.rm=T)+count.sim[j,5]*mean(NCSLdata$distance2, na.rm=T)+count.sim[j,6]*mean(NCSLdata$percentblack, na.rm=T)+count.sim[j,7]*0.75+count.sim[j,8]*mean(NCSLdata$turnout, na.rm=T)+count.sim[j,9]*0.75*med[i]	
		
		}
		
	pe[i]<-mean(invlogit(pc1))
	lb[i]<-quantile(invlogit(pc1), 0.025)
	ub[i]<-quantile(invlogit(pc1), 0.975)
	pe2[i]<-mean(invlogit(pc2))
	lb2[i]<-quantile(invlogit(pc2), 0.025)
	ub2[i]<-quantile(invlogit(pc2), 0.975)	
	
	}

#pdf("PreProbNCD.pdf")	
plot(med, pe, type="l", lty=1, lwd=4, ylim=c(0,1), main="Probability of Being Redistricted From Democrats to Republicans", ylab="Predicted Probability", xlab="Median Household Income", sub="North Carolina Model")
lines(med, lb, lty=2, lwd=2)
lines(med, ub, lty=2, lwd=2)
lines(med, pe2, lty=1, lwd=4, col="darkgrey")
lines(med, lb2, lty=2, lwd=2, col="darkgrey")
lines(med, ub2, lty=2, lwd=2, col="darkgrey")
legend(100, 1.0, c("Incumbent with 50% of Vote", "Incumbent with 75% of Vote", "95% Confidence Intervals"), lty=c(1,1,2), col=c("black", "darkgrey", "black"))
#dev.off()

library(mvtnorm)
library(arm)
sims<-500
count.sim<-rmvnorm(sims, model5.NCr$coef, vcov(model5.NCr))
calc<-400
med<-seq(0, 200, length=calc)
pe<-numeric(calc)
lb<-numeric(calc)
ub<-numeric(calc)
pe2<-numeric(calc)
lb2<-numeric(calc)
ub2<-numeric(calc)
for(i in 1:calc){
	pc1<-numeric(sims)
	pc2<-numeric(sims)
	for(j in 1:sims){
	pc1[j]<-count.sim[j,1]+count.sim[j,2]*med[i]+count.sim[j,3]*mean(NCSLdata$NCHouseDem, na.rm=T)+count.sim[j,4]*mean(NCSLdata$NCHouseRep, na.rm=T)+count.sim[j,5]*mean(NCSLdata$distance2, na.rm=T)+count.sim[j,6]*mean(NCSLdata$percentblack, na.rm=T)+count.sim[j,7]*0.5+count.sim[j,8]*mean(NCSLdata$turnout, na.rm=T)+count.sim[j,9]*0.5*med[i]
	pc2[j]<-count.sim[j,1]+count.sim[j,2]*med[i]+count.sim[j,3]*mean(NCSLdata$NCHouseDem, na.rm=T)+count.sim[j,4]*mean(NCSLdata$NCHouseRep, na.rm=T)+count.sim[j,5]*mean(NCSLdata$distance2, na.rm=T)+count.sim[j,6]*mean(NCSLdata$percentblack, na.rm=T)+count.sim[j,7]*0.75+count.sim[j,8]*mean(NCSLdata$turnout, na.rm=T)+count.sim[j,9]*0.75*med[i]	
		
		}
		
	pe[i]<-mean(invlogit(pc1))
	lb[i]<-quantile(invlogit(pc1), 0.025)
	ub[i]<-quantile(invlogit(pc1), 0.975)
	pe2[i]<-mean(invlogit(pc2))
	lb2[i]<-quantile(invlogit(pc2), 0.025)
	ub2[i]<-quantile(invlogit(pc2), 0.975)	
	
	}
#pdf("PreProbNCR.pdf")	
plot(med, pe, type="l", lty=1, lwd=4, ylim=c(0,1), main="Probability of Being Redistricted From Republicans to Democrats", ylab="Predicted Probability", xlab="Median Household Income", sub="North Carolina Model")
lines(med, lb, lty=2, lwd=2)
lines(med, ub, lty=2, lwd=2)
lines(med, pe2, lty=1, lwd=4, col="darkgrey")
lines(med, lb2, lty=2, lwd=2, col="darkgrey")
lines(med, ub2, lty=2, lwd=2, col="darkgrey")
legend(100, 0.2, c("Incumbent with 50% of Vote", "Incumbent with 75% of Vote", "95% Confidence Intervals"), lty=c(1,1,2), col=c("black", "darkgrey", "black"))
#dev.off()

##########PredProb
library(mvtnorm)
library(arm)
sims<-500
count.sim<-rmvnorm(sims, model5.CA$coef, vcov(model5.CA))
calc<-400
med<-seq(0, 200, length=calc)
pe<-numeric(calc)
lb<-numeric(calc)
ub<-numeric(calc)
pe2<-numeric(calc)
lb2<-numeric(calc)
ub2<-numeric(calc)
for(i in 1:calc){
	pc1<-numeric(sims)
	pc2<-numeric(sims)
	for(j in 1:sims){
	pc1[j]<-count.sim[j,1]+count.sim[j,2]*med[i]+count.sim[j,3]*mean(CASLdata$assdem, na.rm=T)+count.sim[j,4]*mean(CASLdata$assrep, na.rm=T)+count.sim[j,5]*mean(CASLdata$near_dist, na.rm=T)+count.sim[j,6]*mean(CASLdata$perAA, na.rm=T)+count.sim[j,7]*0.5+count.sim[j,8]*mean(CASLdata$turnout, na.rm=T)+count.sim[j,9]*0.5*med[i]
	pc2[j]<-count.sim[j,1]+count.sim[j,2]*med[i]+count.sim[j,3]*mean(CASLdata$assdem, na.rm=T)+count.sim[j,4]*mean(CASLdata$assrep, na.rm=T)+count.sim[j,5]*mean(CASLdata$near_dist, na.rm=T)+count.sim[j,6]*mean(CASLdata$perAA, na.rm=T)+count.sim[j,7]*0.75+count.sim[j,8]*mean(CASLdata$turnout, na.rm=T)+count.sim[j,9]*0.75*med[i]	
		
		}
		
	pe[i]<-mean(invlogit(pc1))
	lb[i]<-quantile(invlogit(pc1), 0.025)
	ub[i]<-quantile(invlogit(pc1), 0.975)
	pe2[i]<-mean(invlogit(pc2))
	lb2[i]<-quantile(invlogit(pc2), 0.025)
	ub2[i]<-quantile(invlogit(pc2), 0.975)	
	
	}

#pdf("PreProbCAD.pdf")	
plot(med, pe, type="l", lty=1, lwd=4, ylim=c(0,0.5), main="Probability of Being Redistricted From Democrats to Republicans", ylab="Predicted Probability", xlab="Median Household Income", sub="California Model")
lines(med, lb, lty=2, lwd=2)
lines(med, ub, lty=2, lwd=2)
lines(med, pe2, lty=1, lwd=4, col="darkgrey")
lines(med, lb2, lty=2, lwd=2, col="darkgrey")
lines(med, ub2, lty=2, lwd=2, col="darkgrey")
legend(100, 0.5, c("Incumbent with 50% of Vote", "Incumbent with 75% of Vote", "95% Confidence Intervals"), lty=c(1,1,2), col=c("black", "darkgrey", "black"))
#dev.off()

library(mvtnorm)
library(arm)
sims<-500
count.sim<-rmvnorm(sims, model5.CAr$coef, vcov(model5.CAr))
calc<-400
med<-seq(0, 200, length=calc)
pe<-numeric(calc)
lb<-numeric(calc)
ub<-numeric(calc)
pe2<-numeric(calc)
lb2<-numeric(calc)
ub2<-numeric(calc)
for(i in 1:calc){
	pc1<-numeric(sims)
	pc2<-numeric(sims)
	for(j in 1:sims){
	pc1[j]<-count.sim[j,1]+count.sim[j,2]*med[i]+count.sim[j,3]*mean(CASLdata$assdem, na.rm=T)+count.sim[j,4]*mean(CASLdata$assrep, na.rm=T)+count.sim[j,5]*mean(CASLdata$near_dist, na.rm=T)+count.sim[j,6]*mean(CASLdata$perAA, na.rm=T)+count.sim[j,7]*0.5+count.sim[j,8]*mean(CASLdata$turnout, na.rm=T)+count.sim[j,9]*0.5*med[i]
	pc2[j]<-count.sim[j,1]+count.sim[j,2]*med[i]+count.sim[j,3]*mean(CASLdata$assdem, na.rm=T)+count.sim[j,4]*mean(CASLdata$assrep, na.rm=T)+count.sim[j,5]*mean(CASLdata$near_dist, na.rm=T)+count.sim[j,6]*mean(CASLdata$perAA, na.rm=T)+count.sim[j,7]*0.75+count.sim[j,8]*mean(CASLdata$turnout, na.rm=T)+count.sim[j,9]*0.75*med[i]	
		
		}
		
	pe[i]<-mean(invlogit(pc1))
	lb[i]<-quantile(invlogit(pc1), 0.025)
	ub[i]<-quantile(invlogit(pc1), 0.975)
	pe2[i]<-mean(invlogit(pc2))
	lb2[i]<-quantile(invlogit(pc2), 0.025)
	ub2[i]<-quantile(invlogit(pc2), 0.975)	
	
	}

#pdf("PreProbCAR.pdf")	
plot(med, pe, type="l", lty=1, lwd=4, ylim=c(0,1), main="Probability of Being Redistricted From Republicans to Democrats", ylab="Predicted Probability", xlab="Median Household Income", sub="California Model")
lines(med, lb, lty=2, lwd=2)
lines(med, ub, lty=2, lwd=2)
lines(med, pe2, lty=1, lwd=4, col="darkgrey")
lines(med, lb2, lty=2, lwd=2, col="darkgrey")
lines(med, ub2, lty=2, lwd=2, col="darkgrey")
legend(0, 1.0, c("Incumbent with 50% of Vote", "Incumbent with 75% of Vote", "95% Confidence Intervals"), lty=c(1,1,2), col=c("black", "darkgrey", "black"))
#dev.off()


